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SUMMARY 


Two relaxation schemes for Chebyshev spectral multigrid methods are presented for elliptic 
equations with Dirichlet boundary conditions. The first scheme is a pointwise-preconditioned 
Richardson relaxation scheme and the second is a line relaxation scheme. The line relaxation 
scheme provides an efficient and relatively simple approach for solving two-dimensional spectral 
equations. Numerical examples and comparisons with other methods are given. 


INTRODUCTION 


For limited-area problems with general (non-periodic) boundary conditions, Chebyshev spectral 
methods give exponential convergence for smooth solutions. However, except in some very simple 
cases (e.g., one-dimensional constant-coefficient problems), Chebyshev approximations usually lead 
to full linear systems which cannot be solved efficiently by direct methods, and iterative methods 
must be used. Unfortunately, designing efficient iterative methods for discrete spectral equations 
has proven difficult, especially for problems with non-constant coefficients (ref. 1). Perhaps the 
most promising technique to date for solving spectral discretizations of elliptic problems is the 
spectral multigrid method (ref. 2, 3). However, the best relaxation schemes known today are 
complicated to apply. In this paper we introduce two simpler relaxation schemes and investigate 
their performance. 


As prototype problems we consider one- and two-dimensional elliptic equations with Dirichlet 
boundary conditions on simple geometric domains. In one dimension we consider 


-u"(x) = f(x), 
u(+l) = a, 

The two-dimensional prototype problem is 

-A u(x,y) = f{x,y), 
u{x,y) = g(x,y), 


\x\ < 1, 
u(— 1) = b. 


M,|y| < i, 

1*1 = i.M = L 


(1) 

(2) 


We discretize these problems by Chebyshev collocation. For example, for the two-dimensional 
problem (2), the solution u(x, y) is approximated by a set of discrete values uj^ on the Chebyshev 
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grid {( Xj,Vk ) = (cos(jn/N x ),cos(kTr/N y )) | 0 < j < N x , 0 < k < N y }, with the requirement that 
problem (2) be satisfied on this grid, i.e., 

+vj! ) ) = *). 

u j,k = 9{ x j 

where and u - >. are values of the second-order derivatives of the Chebyshev approximation 

Y,m=o ^n=0 UmnTm(x)T n (y ) to u(x, y) on the Chebyshev grid. For simplicity, we will assume here 
that N x — N y = N] however, the codes described in this paper do not require this. 

The discrete problem (3) can be expressed in form of a linear system 


1 < j < N x , 0 < k < N y 
j = 0,j = N x ,k =0,k = N y 


AU = F 


( 4 ) 


Unfortunately, the matrix A, formulated by Chebyshev collocation approximations, is full and 
non-symmetric. For two-dimensional problems, direct methods (like Gaussian elimination) would 
require 0(N G ) operations for factorization and 0(iV 4 ) for the subsequent solution, which is far too 
much work to be practical. Thus, iterative methods must be used. 

1 

I 


THE POINTWISE PRECONDITIONED RICHARDSON RELAXATION SCHEME 

The most efficient method available today for solving (4) and its generalizations to other 
elliptic problems is the spectral multigrid method of Zang et al. (ref. 2, 3), which employs finite- 
difference preconditioned Richardson iteration as the relaxation scheme in a multigrid context. 
Preconditioned Richardson relaxation for (4) takes the form 

V i — V + uH(F - AV), (5) 

where V is the current approximation to U,u is a relaxation parameter, and H is the 
preconditioner. The criteria for choosing a preconditioner H are: 

• H should give fast multigrid convergence, 

• H should be easy and cheap to generate or apply. 

The finite-difference preconditioning of Zang et al. (ref. 2, 3) gives fast convergence, but applying it 
requires solving (or nearly solving) a finite-difference discretization on the nonuniform Chebyshev 
grid. This procedure is complicated and expensive. Are there alternatives which are simpler 
and still effective? Achi Brandt (personal communication, 1983) has suggested that pointwise 
preconditioning based on the (variable) Chebyshev mesh spacing might work well. In this section, 
we investigate the performance of this simple preconditioner when applied to the problem (4). 
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IllilUIIIIIUli 


The One-Dimensional Case 


Formulation 

As an analogue of the Gauss-Seidel relaxation for a finite-difference method, the pointwise 
preconditioning for the Chebyshev discretization takes the form 

hj / fi x 

Vj « — Vj+u-±rj, (b) 

where hj — ( x j — 1 — xj-^i)/2 is the effective grid size at the point x j , tj is the the residual 
R = F — AV at Xj, and u is a relaxation parameter to be chosen to accelerate the convergence. 
Note that (6) is equivalent to choosing the preconditioning matrix H in (5) as a diagonal matrix 

= (7) 


Analysis 

The evolution of the error E = V - U in the Richardson relaxation (5) is described by 

E i — (7 - u>HA)E. 


( 8 ) 


Therefore, the convergence factor for (5) on a single grid is 

&SG — p{I ~ wHA), 

where p denotes the spectral radius. Likewise, the multigrid smoothing factor for (5), when used as 
a smoother in a multigrid method (e.g., ref, 4), is 

Jl = p(G{I-u;HA)), (9) 

where G represents the perfect coarse-grid correction, i.e., set all low modes of the error to zero. 


For the simple preconditioning (7), our numerical computations show that the eigenvalues of the 
matrix HA are all positive real numbers. The maximum eigenvalue is A max « 5.0, the middle is 
A m i<j ~ 1.5, and the minimum is A m ; n « 0(N~ 2 ). The formulas of Zang et. al. (ref. 2, 3) then give a 
good approximation to the optimal u> and p , , namely, 


u 


T -^mid 


0.325, 




■^mid 


T Amid 


0 . 6 . 


( 10 ) 


Indeed, computing the smoothing factor directly from (9) using u> — 0.325, we find that p < 0.6 for 
all N < 512. 


To take into account the effects of grid transfers (omitted in the smoothing analysis above), 
we use the following two-grid analysis. The evolution of the error E in one two-grid V (ni, re- 
cycle (where ni and ri 2 specify the number of relaxation sweeps before and after the coarse-grid 
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correction, respectively) is described by the matrix 

T = (I - uHA f ) n2 {I - PA~ l RA)(I - uHA f ) n K (11) 

Here, R represents the fine-to-coarse grid transfer (we use injection), P represents the coarse-to- 
fine grid transfer (we use Chebyshev interpolation), and Af and A c represent the discrete operator 
matrix in (4) on the fine and coarse grids, respectively. Note that (11) assumes that the coarse-grid 
problem is solved exactly. 

We computed the two-grid convergence factor <jtg = P (T) for N < 512 using different 
values of ui, and the numerical results show that u> = 0.325 again gives the optimal convergence 
factor (or very close to it). Using that constant value, we find that the smoothing factor per sweep 
Ps = (<7TG) 1 ^ ni+n2 ^ satisfies 

0.5 < fx 3 < 0.6 

for all N < 512. A similar analysis for the one-dimensional Helmholtz problem 

A u(x) - u"{ x) = f(x) (12) 

shows that with various choices of A and boundary conditions (Dirichlet, Neumann and mixed), an 
appropriate pointwise preconditioner also yields the smoothing factor per sweep < 0.6. 

We have developed FORTRAN-77 routines to implement the Chebyshev multigrid method 
using the pointwise preconditioner as described above. The code has been used to solve the 
problem (12) with various choices of u(x), A, and boundary conditions. The observed convergence 
factor per sweep /j, 3 is smaller than 0.60 for all cases tested, in agreement with the analysis 
presented above. 


The Two-Dimensional Case 


Formulation 


We note that Gauss-Seidel relaxation for the second-order centered finite difference 
approximation to (2) can be written as 



Uj, k + 


h 2 

4 r iA> 


where rjj. is the finite-difference residual. A natural analogue for the Chebyshev collocation 
discretization (3) is 


Vj.k < — Vj,k + v 


.2 !h) + 2 1 hi 


r j,k> 


(13) 


where hj and are the grid sizes at the point (xj, y*.), fj k := fj k — [— (uj z fc z ^ + Vj V ^)] is the 
residual of Chebyshev discretization, and oj is a relaxation parameter to be chosen to accelerate 
the convergence. Clearly, the iteration (13) is a special case of the Richardson iteration (5), with 
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a diagonal preconditioner H with diagonal entries {H)jkjk — + 2//i|) 1 . This preconditioner 

is easy and fast to apply. Does it gives a fast convergence? Unfortunately, the following analysis 
shows that the answer is no. 


Analysis 

Computational results indicate that the eigenvalues of the matrix HA are all positive real 
numbers. Again, good approximations to the optimal u and /i can be obtained by 


u 


^max T Aq Ua 




Wax ^qua 


+ A 


qua 


(14) 


where A max is the maximum eigenvalue and A qua is the quarter eigenvalue (ref. 1). More precise 
values of the optimal u> and /z, can be obtained by actually computing the spectral radius 
p(G(I - wH A)) for different choices of u> and comparing the results. For N < 32, the eigenvalues 
Amax and A qua , u and p computed by (14) and the optimal u; and p are listed in Table I. Since p is 
large and increases with iV, these results suggest that the pointwise preconditioner (13) will not be 
a good multigrid smoother. 

Table I also lists the two-grid smoothing factors per sweep fi s = (p(T)) 1 / (ni+ ” 2) computed from 
the matrices in (11) for N < 32 using u = 0.36. These results again show that the pointwise 
preconditioning (13) does not give fast convergence. 

We have implemented the pointwise preconditioning (13) in a multigrid solver written in 
Fortran 77. Computational results from a number of test cases confirm the above analysis: we 
conclude that the pointwise preconditioning does not give fast convergence. 


Table I. Multigrid Analysis of Two-Dimensional Pointwise Preconditioning 


N 

Eigenvalues of HA 

By (14) 

By computation 

Amax 

Aqua 

u 

At 

^ opt 

A 1 

At, 

4 

3.00 

1.83 

0.41 

0.24 

0.35 

0.28 

0.51 

8 

4.10 

1.26 

0.37 

0.53 

0.35 

0.52 

0.68 

16 

4.57 

0.95 

0.36 

0.66 

0.36 

0.75 

0.80 

32 

4.76 

0.78 

0.36 

0.72 

0.36 

0.82 

0.88 
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THE LINE RELAXATION METHOD 


The poor performance of pointwise preconditioning in two dimensions can be understood in 
terms of the anisotropy introduced by the nonuniform Chebyshev collocation grid. Since the mesh 
spacing varies with x and y, at any given point (x, y ) the coupling in the discrete operator in (3) 
may be stronger in x or in y. In finite-difference multigrid methods, point relaxation performs ] 

poorly in such anisotropic cases, and the cure is to use alternating direction line relaxation. Thus, 
it is reasonable to try an analogous approach for the Chebyshev discretization. 


Formulation 


To formulate the line relaxation method, we express the discrete problem (3) in the matrix form 


(n + V)U=F, 


(15) 


where W, and V correspond to the horizontal part (— d 2 jdx 2 ) and vertical part (—d^/dy 2 ) of the 
Laplacian operator, respectively. Starting from an approximation K old to the solution U , one sweep 
of (alternating direction) line relaxation based on (15) consists of the following two parts: 


1. Sweep along the x-direction. On each grid line parallel to x-axis, use the values of V old except 
those on the current line, and solve for values on the current line by solving (15). This can be 
expressed in the matrix form as 


(H + V d )F mid — F — V Q V old , (16) 

where Vd and V 0 denote the diagonal and off-diagonal parts of the matrix V, respectively. Note 
that the entries of are known (ref. 1) and Vd is a constant on each grid line parallel to the x- 
axis. Thus, the system (16) can be decoupled into (A T — 1) one-dimensional discrete problems, 
each of which is a Chebyshev collocation approximation to a Helmholtz equation on an interior 
grid line parallel to x-axis; the x-directional sweep consists of solving these equations. 

2. Sweep along the y-direction. The y-direction sweep is basically the same as the x-direction 
sweep except that we now work on grid lines that are parallel to y-axis and use values of V mid 
instead of V old . The equation we need to solve is 


(H d + V)V"' :W = F- 


( 17 ) 


where 77j and H 0 are the diagonal and off-diagonal parts of H. As in the x-direction sweep, the 
two-dimensional problem (17) is solved by solving (N - 1) one-dimensional Helmholtz equations. 

It turns out that as it stands, the line relaxation (16)-(17) is not a good multigrid smoother; 
however, this can be fixed as follows. Let C mid = V mid — V°* d and C new = I/ new — V rmid denote the 
corrections for V old and l /mid , and R old = F — AV' M and R mid = F — AV mid denote the residuals 
of V old and V mid , respectively. Rewriting equations (16) and (17) as correction equations and 
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introducing a relaxation parameter u> (to be determined by analysis to accelerate the convergence), 
we obtain 

(H + V d )C mid = uR o]d , (U d + V)C new = b >R mid (18) 

We refer to (18) as the collocation version of the line relaxation method. 

It is not practical to implement the collocation version because there are no fast solvers available 
for the collocation approximations, even for one-dimensional problems. However, in the multigrid 
context, a relaxation scheme functions as a smoother rather than a solver: instead of solving 
each problem exactly, we only need to smooth out the error, i.e., reduce high modes in the error. 
Therefore, it is reasonable to replace the one-dimensional problems in (18) by approximate versions 
which can be solved efficiently. We consider two alternatives as follows. 

In the first, we replace the collocation discretizations of the one-dimensional Helmholtz 
equations in (18) by tau discretizations. Tau approximations have the same exponential 
convergence as collocation method, but can be solved directly in 0(7Vlog N) operations. This leads 
to the tau version of the line relaxation method, and the total work of one x or y-direction sweep is 
0(N 2 log N ). As we will see below, this tau version turns out to be an efficient multigrid smoother. 

In the second, we replace the collocation discretizations of the one-dimensional Helmholtz 
equations in (18) by finite-difference discretizations. This leads to the finite-difference version of 
the line relaxation method, which has two obvious advantages over the tau version. First, it is 
faster because it eliminates the transforms required in tau version, thus reducing the operation 
count for solving each one-dimensional problem from 0(N\ogN) to O(N). Second, it can be 
extended to solve more generalized problems, e g., problems with variable coefficients. As we will 
see below, this finite-difference version also turns out to be an efficient multigrid smoother, even in 
the case of variable coefficients. 


Analysis 


As in the case of the pointwise preconditioned Richardson relaxation, we can analyze the 
performance of the line relaxation methods described above by computing the eigenvalues of the 
corresponding interation matrices. Because the tau version cannot be expressed in matrix form 
like (18), we will only do the analysis for the collocation and finite-difference versions. Note that 
the tau and collocation versions are nearly the same, so the analysis for collocation version should 
give a good prediction for the performance of the tau version. In this section, we will give details of 
the analysis for finite-difference version and only list results for collocation version. 


Smoothing Analysis 

For the finite-difference version of the line relaxation iteration, the error evolution is described 

b y 


£mid ^ [/ _ 

• w{u fd + v d )-\n + V)]E° U , 

(19) 

£new < [7 _ 

u>{U d + V fd )~ x {U + V)]£ mid . 

(20) 
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where 'H fd and V fd are the finite-difference analogues of the collocation discretization matrices U 
and V, respectively. Therefore, the error evolution matrix for one relaxation is 

S = [I- u(H d + V fd )~\U + V)][/ - a ,(U fd + V d )- l {U 4- V)]. (21) 

The matrices = {U fd + V d ) _1 0H + V) and S v = {U d + V fd )- l {U + V) have the same 
eigenvalues (since x and y can be interchanged in the Laplacian operator), so we can focus on 
just the x-direction sweep (19). The eigenvalues of S% are all positive real numbers, so we can 
use formulas (14) to obtain approximate values of lj and Jl (squaring p to represent the effect of 
both the x and y sweeps). These values are listed in Table II for N < 32, along with the optimal 
relaxation parameter u) and corresponding multigrid smoothing factor p = p{GS) computed 
directly. These results suggest that for large values of truncation number N , u opt « 0.6 and 
jX < 0.5, independent of the grid size. Corresponding results for the collocation version are listed 
in Table III. 


Multigrid Analysis 

For a multigrid V{n\, ri 2 )-cycle, if we use zeros as initial guesses on all coarse grids (which is 
a natural choice because the coarse-grid solution is a correction to the solution on the next finer 
grid), then we can write out the error evolution matrix explicitly as 

M = S n2 [I - PGR(H + V)]5 ni . (22) 

This represents a procedure of ni pre-relaxations ( S ni ) followed by a coarse-grid-correction 
(/ _ PGR (W. + V)) and then n 2 post-relaxations (S n2 ). The matrix S is the error evolution 
matrix of one relaxation on the finest grid defined in (21). The central part I — PGR{H + V) 
represents the coarse-grid-correction, where R represents the fine-to-coarse grid transfer (we use 
injection) and P represents the coarse-to-fine grid transfer (we use Chebyshev interpolation). The 
matrix G is defined on the next coarser grid as follows: on the coarsest grid, G = (H + V) -1 (which 
means the coarsest grid problem is solved exactly); otherwise, 

G = [I - M] * + V) -1 , (23) 

which represents a multigrid solution procedure on that grid. Note that (23) is actually a recursive 
definition, since the matrix M in (23) includes another matrix G on the next coarser grid. 

Tables II and III also list computed values of smoothing factor per sweep p s - (p(M)) 1 ^ ni+n ^ 
for the case lo = 0.6, ni = 2, and n 2 — 1. These results suggest that the smoothing factor of the 
line relaxation method is less than 0.5, independent of the grid size. Note that while we could also 
use Chebyshev restriction instead of injection for the fine-to-coarse grid transfer R , our numerical 
experience shows very little difference between these two choices. 
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Table II. Analysis of the Finite-Difference Version 


N 

Eigenvalues of 5-// 

By (14) 

By computation 

^max 

^qua 

u 

V 

^opt 

V 

Vs 

4 

1.995 

1.000 

0.669 

0.110 

0.58 

0.110 

0.181 

8 

2.513 

1.000 

0.569 

0.186 

0.60 

0.168 

0.293 

16 

2.780 

0.995 

0.530 

0.224 

0.60 

0.271 

0.364 

32 

2.898 

0.815 

0.539 

0.315 

0.60 

0.366 

0.421 


Table III. Analysis of the Collocation Version 


N 

Eigenvalues of Sy, 

By (14) 

By computation 

^max 

^qua 

UJ 

V 

t^opt 

V 

Vs 

4 

1.651 

1.000 

0.754 

0.060 

0.68 

0.120 

0.302 

8 

2.322 

0.922 

0.616 

0.186 

0.60 

0.216 

0.328 

16 

2.701 

0.810 

0.570 

0.290 

0.58 

0.326 

0.380 

32 

2.869 

0.700 

0.560 

0.370 

0.60 

0.410 

0.428 


Computational Results 


We have implemented the tau and finite-difference versions of the line relaxation scheme 
described above in a Chebyshev collocation multigrid solver for the two-dimensional Helmholtz 
problem 

Xu(x, y) - A u(x, y ) = f(x, y), \x\, |y| < 1, 

u(x,y) = g(x,y), 1x1 = 1 , 12/1 = 1, 

with various choices of /, g, and A. For both versions, the observed convergence factor per sweep is 
less than 0.5 for all cases tested, in agreement with the analysis above. The finite-difference version 
turns out to have slightly better convergence factors than the tau version, but the difference is 
minor. 
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Comparisons with Other Methods 


In this section we compare the line relaxation spectral multigrid method developed above to two 
other methods for solving the two-dimensional prototype problem (2). The first is a conventional 
finite-difference multigrid method; the second is a matrix diagonalization technique. We do not 
compare with the method of Zang et. al. (ref. 3) since the details presented in that paper were not 
enough to allow programing the method. All computations are done on a SUN SPARCstation2 
using double precision; the machine round-off error is about 2.22 x 10 


Conventional Finite-Difference Multigrid Method 

The finite-difference discretization is the usual second-order five-point scheme on a uniform 
grid. The finite-difference multigrid method uses Gauss-Seidel (Red-Black) iteration as a relaxation 
scheme, the fine-to-coarse grid transfer is half-injection, the coarse-to-fine grid transfer is bilinear 
interpolation, and the multigrid V-cycle algorithm is used. 

According to computations, the average execution time of one U(2, l)-cycle of the finite- 
difference multigrid method is approximately (0.56 x 1CT 4 ) N 2 seconds, and (0.21 x 10~ 3 ) N 2 log 2 N 
seconds for line relaxation spectral multigrid method. Therefore, for the same grid sizes, one 
V(2, l)-cycle of the finite-difference multigrid method is approximately 3.75 log 2 N times faster 
than the line relaxation spectral multigrid method. 

However, because spectral methods have exponential convergence and finite-difference 
methods only have polynomial convergence, when high accuracy is required, finite-difference 
multigrid methods must use much bigger grid sizes than spectral methods. The result is that 
the line relaxation spectral multigrid method is faster than finite-difference when high accuracy 
is required. As a specific example, consider the prototype problem (2) with true solution 
u(x, y) - e 2x+y cos(n(x + 4 y + 0.25)). The relation between accuracy and execution time required 
to achieve that accuracy is plotted in Figure 1 for both methods. We can see that when low 
accuracy is required, the finite-difference multigrid method is much faster than the line relaxation 
spectral multigrid method, but the situation is reversed when high accuracy is required. The 
crossover point for this problem is at an accuracy of about one percent error. The same conclusion 
would hold for finite-difference methods of higher (fixed) orders, although the crossover point 
would shift. Variable-order finite-difference methods could be expected to perform more like the 
spectral method, at a cost of considerable complexity. 


Matrix Diagonalization Technique 

The matrix diagonalization technique is introduced in (ref. 5) as a direct solver for the 
Chebyshev spectral approximation to the Poisson equation with Dirichlet boundary conditions. 
This technique requires a preprocessing step, which involves computing the eigenvalues and 
eigenvectors of a one-dimensional operator matrix (0(A 3 ) operations), and a solution step, which 
involves one-dimensional matrix multiplications ( 0(N 3 ) operations). 
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Figure 1. Execution time: LR-SMG vs FD-MG 


To compare execution times, we note that the line relaxation spectral multigrid method 
usually takes approximately 10 F-cycles to solve to the level of machine precision. Thus, Figure 2 
compares the execution time of 10 V-cycles of the line relaxation spectral multigrid method with 
the execution time of solving the same problem directly by the matrix diagonalization method 
(including the preprocessing step). These results show that the matrix diagonalization method is 
quite fast for small grid sizes, but as the grid size grows, it becomes slower than the line relaxation 
spectral multigrid method. This is because the line relaxation spectral multigrid method is an 
0(N 2 log N) method, while the matrix diagonalization method requires 0(N 3 ) operations (even 
without the preprocessing step). 

The matrix diagonalization technique is very efficient for problems with constant coefficients, 
especially when repeated solutions are required. However, this technique can only handle problems 
with constant coefficients. As shown below, the line relaxation spectral multigrid method is able to 
solve problems with non-constant coefficients. 
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Figure 2. Execution time: LR-SMG vs MD 


Extension to Problems With Variable Coefficients 


As a test problem with variable coefficients we consider 

~M? {X ’ v) lL' ,{X ' y) ) - Ty = h,v), 

u(x,y) = g{x,y), 

where the coefficient functions and the true solution are 

a(x, y) — b(x, y) — 1 + £e cos ^ n ^ x+y ^\ 


u(x, y) = sin(a7rx + sin(a7ry 4- j). 


W»|y| < i» 

1*1 = i Av\ = i. 


(24) 


(25) 


( 26 ) 


The parameter e measures how far the coefficients are away from the constant 1, (3 measures the 
oscillation of the coefficients, and a measures the oscillation of the solution. 
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Implementation of the Line Relaxation Spectral Multigrid Method 

The implementation of the finite-difference version of the line relaxation method is basically the 
same as for the constant coefficient case except for the following: 


1. On each grid line, the one-dimensional problem is not a Helmholtz equation anymore. For 
example, on a gird line y — y k which is parallel to x-axis, we now solve a problem like 

~lL " 11 * V d{x,yk)v(x,y k ) = h(x,y k ) (27) 

by using a. second-order finite-difference approximation on the Chebyshev grid. 

2. To compute values of note that the interior equation in (24) can be rewiitten as 


d 2 u dadu d 2 u 

a dx 2 dx dx dy 2 


db du 
dy dy 


x\,\y\ < 1 


and the Chebyshev collocation approximation to (28) can be written as 

{-AV XX - A X V X } U - {-BVyy - ByVy } U=F , 


(28) 


(29) 


where A and B axe diagonal matrices containing the values of the coefficients a(xj, y k ) and 
b(xj,y k ), A x and B x are diagonal matrices containing the values of the derivatives ^a(xj,y k ) 
and Jf-b(xj,y k ) (which can be computed from values a(xj, y k ) and b(xj,y k )), and V x , V xx , 

V y , and V yy are Chebyshev differentiation matrices. Therefore, % — —AD XX — A X V X and 
V = -BVyy - ByVy ; generating the diagonal entries of U and V is straightforward. 

3. On coarse grids, we need to use so-called “filtered” coefficients a(x, y) and b(x, y) to formulate 
the coarse grid problems; i.e., the coefficients a(x,y ) and b(x,y ) are evaluated on the finest grid 
and then transferred to the coarser grids by Chebyshev restriction (ref. 3). 


Computational Results 

We have run the line relaxation spectral multigrid method for different values of parameters e, a 
and p. For a = 1.0 and N x = N y = 32, the smoothing factor is graphed in Figure 3 as a function of 
e and a. Here we have chosen to measure the smoothing by the “smoothing factor per work unit” 
defined by fi w — (r 2 /Vi) ro//r , where r\ and r 2 are residual norms before and after one multigrid 
V-cycle, r is the execution time of one cycle and r 0 is the execution time of one relaxation. These 
results show that for a wide range of e and /?, the method converges relatively quickly. 

I 11 (ref. 3) the same test problem (24) was solved using the Richardson relaxation (5) using 
two-dimensional finite- difference preconditioning; incomplete LU decomposition was used to 
approximately solve the finite difference approximation on the Chebyshev grid. With only limited 

details of the formulation and results of this method, it is difficult to make a complete comparison 

to the line relaxation method considered here. However, it appears that the line relaxation method 
gives convergence factors at least as small as those in (ref. 3); moreover, it is simpler. 
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CONCLUSIONS 


The pointwise preconditioning is simple and fast to apply. It is very efficient for one-dimensional 
problems. Unfortunately, it does not give fast multigrid convergence for two-dimensional problems. 

The line relaxation method provides a new approach to accelerate the multigrid Chebyshev 
spectral method for solving two-dimensional elliptic problems. It is efficient (yielding multigrid 
smoothing factors no larger than 0.5 per sweep) and inexpensive (requiring 0{N 2 log .V) operations 
per sweep). . •• - 

When high accuracy is required, the spectral multigrid method using line relaxation is orders 
of magnitude faster than a conventional finite-difference multigrid method, due primarily to the 
exponential convergence of the spectral discretization. Compared to other methods for solving 
the discrete spectral equations, the line relaxation method also has advantages: it is comparable 
in efficiency to matrix diagonalization and finite-difference preconditioned Richardson relaxation, 
but can solve problems with variable coefficients which the former cannot, and is simpler than the 
latter. 
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